PCLE Replication
Emorie Beck
04/04/2018
Workspace
Packages
library(survey)
library(mi)
library(psych)
library(MatchIt)
library(lavaan)
library(semTools)
library(rstan)
library(lme4)
library(brms)
library(rstanarm)
library(tidybayes)
library(MuMIn)
library(parallel)
library(gridExtra)
library(knitr)
library(kableExtra)
library(stargazer)
library(plyr)
library(stringr)
library(haven)
library(tidyverse)
library(ggridges)
data_path <- "~/Box/network/other projects/PCLE Replication"
data_path <- "https://github.com/emoriebeck/life_events/raw/master"
model_path <- "~/Box/Models/PCLE Replication"Data
Load Raw Data
First we’ll load in the Codebook that lists which dataset each variable is in, what’s it’s named in those data files, what I’d like to name it for consistency, the scale it’s measured on, the year it was measured, whether it was reverse keyed, whether I want to recode the variable, and what rules I want to use to create composites (for matching variables).
meta <- readxl::read_xlsx(sprintf("%s/data/Codebook.xlsx", data_path)) %>%
mutate(Item = stringr::str_to_lower(Item))
all.old.cols <- (meta %>% filter(class == "proc" & Year == 0))$Item
all.new.cols <- (meta %>% filter(class == "proc" & Year == "0"))$new_name
# create short function to read in separate files for each wave
read_fun <- function(file, year){
print(year)
old.names <- (meta %>% filter(Year == year & class %in% c("group", "predictor", "proc")))$Item
new.names <- (meta %>% filter(Year == year & class %in% c("group", "predictor", "proc")))$new_name
z <- haven::read_sav(sprintf("%s/data/sav_files/%sp.sav", data_path, file)) %>%
left_join(haven::read_sav(sprintf("%s/data/sav_files/%spkal.sav", data_path, file))) %>%
left_join(haven::read_sav(sprintf("%s/data/sav_files/%sh.sav", data_path, file))) %>%
select(one_of(all.old.cols), one_of(old.names)) %>%
setNames(c(all.new.cols, new.names)) %>%
mutate_all(funs(mapvalues(., seq(-1,-7,-1), c(NA, 0, rep(NA,5)), warn_missing = F))) %>%
group_by(PROC_SID) %>%
mutate(LE_ParDied = max(LE_MomDied,LE_DadDied, na.rm = T),
LE_ParDied = ifelse(is.nan(LE_ParDied) == T, NA, LE_ParDied)) %>%
gather(key = new_name, value = value, -PROC_SID, -PROC_household, -Dem_DOB, -Dem_Sex) %>%
left_join(meta %>% filter(Year == year) %>% select(new_name, rev_code)) %>%
mutate(value = ifelse(rev_code == 0 | is.na(rev_code), value,
reverse.code(keys=-1, items=value, mini=1, maxi=8))) %>%
select(-rev_code)
}
dat <- tibble(
Year = as.character(seq(2005, 2015,1)),
file = c(letters[22:26], paste("b", letters[1:6], sep = ""))) %>%
mutate(data = map2(file, Year, read_fun)) %>%
unnest(data) %>%
group_by(PROC_SID) %>%
mutate(
Dem_DOB = max(Dem_DOB, na.rm = T),
Dem_DOB = ifelse(is.infinite(Dem_DOB) == T, NA, Dem_DOB),
Dem_Sex = max(Dem_Sex, na.rm = T),
Dem_Sex = ifelse(is.infinite(Dem_Sex) == T, NA, Dem_Sex)
)load(url(sprintf("%s/results/data.RData", data_path)))Clean BFI Data
Let’s clean the personality data. to do this, we’ll want to recode the years into waves, figure out what the first wave we have for each person is, and recode the sex and age variables for later. Then, I’ll do some transformations of the data into wide format for later and for matching.
bfi_dat <- dat %>% ungroup() %>%
separate(new_name, c("type", "Item"), sep = "_") %>%
filter(type == "BF") %>%
mutate(wave = as.numeric(mapvalues(Year, seq(2005,2013,4), 1:3))) %>%
group_by(PROC_SID) %>%
mutate(fy = min(Year, na.rm = T),
wave = ifelse(fy != 2005, wave-min(wave) + 1, wave)) %>%
# mutate(wave = seq(1, n(), 1)) %>%
# find people who didn't do any Big 5 responses
group_by(PROC_SID, Item) %>%
mutate(na = sum(!is.na(value))) %>%
ungroup() %>%
# recode gender & center age at first BFI wave (2005)
mutate(sex12 = mapvalues(Dem_Sex, c(1,2), c(1,0), warn_missing = F),
sex.c = as.numeric(scale(sex12, center = T, scale = F)),
age = as.numeric(fy)-Dem_DOB,
age.c = age - mean(age, na.rm = T),
age.c2 = age.c^2, #agec2 = agec^2,
age.c3 = age.c^3) #agec3 = agec^3,
bfi_wide <- bfi_dat %>%
filter(na > 1) %>%
separate(Item, c("Trait", "Item"), 1) %>%
unite(Item, wave, Item, sep = "_") %>%
mutate(Item = sprintf("T%s", Item)) %>%
select(PROC_SID, Trait, Item, value, sex12:age.c3) %>%
spread(key = Item, value = value)
bfi_match <- bfi_dat %>%
filter(na > 1) %>%
unite(Item, Item, wave, sep = "_") %>%
select(PROC_SID, Item, value) %>%
spread(key = Item, value = value)BFI Scale Reliability
We want to make sure our scales are performing as expected, so we’ll calculate Cronbach’s alpha for each scale for each year.
alpha_fun <- function(df){
df <- df %>% select(-PROC_SID)
psych::alpha(df)$total$raw_alpha
}
bfi_wide %>% select(PROC_SID, Trait, T1_1:T3_3) %>%
gather(key = item, value = value, T1_1:T3_3) %>%
separate(item, c("wave", "item"), sep = "_") %>%
spread(key = item, value = value) %>%
group_by(wave, Trait) %>%
nest() %>%
mutate(alpha = map(data, alpha_fun)) %>%
unnest(alpha, .drop = T) %>%
spread(key = wave, value = alpha) %>%
kable(., "html", digits = 2, escape = F, booktabs = T,
col.names = c("Trait", "Wave 1", "Wave 2", "Wave 3"),
caption = "Cronbach's Alpha Scale Reliabilities for the BFI-S Subscales") %>%
kable_styling(full_width = F)| Trait | Wave 1 | Wave 2 | Wave 3 |
|---|---|---|---|
| A | 0.51 | 0.50 | 0.49 |
| C | 0.62 | 0.59 | 0.57 |
| E | 0.65 | 0.65 | 0.66 |
| N | 0.61 | 0.62 | 0.62 |
| O | 0.62 | 0.62 | 0.62 |
Clean Life Event Data
For our life event groups, we need to figure who experienced what. We care less about when as long as it was between 2006 and 2015 and didn’t occur in 2005. Some of the original variables are on different scales, so we’ll recode those. This will eventually leave us with unique codes for each of the life events, where each participant will be coded as “1” if they experienced that event and “0” otherwise.
event_fun <- function(df, event){
print(event)
print(unique(df$Year))
z <- df %>%
select(-type, -file) %>%
spread(key = Year, value = value) #%>%
# mutate(
# `2006` = ifelse(`2005` != 1 & `2006` == 1, 1, 0),
# `2007` = ifelse(`2006` != 1 & `2007` == 1, 1, 0),
# `2008` = ifelse(`2007` != 1 & `2008` == 1, 1, 0),
# `2009` = ifelse(`2008` != 1 & `2009` == 1, 1, 0),
# `2010` = ifelse(`2009` != 1 & `2010` == 1, 1, 0),
# `2011` = ifelse(`2010` != 1 & `2011` == 1, 1, 0),
# `2012` = ifelse(`2011` != 1 & `2012` == 1, 1, 0),
# `2013` = ifelse(`2012` != 1 & `2013` == 1, 1, 0))
# if(event != "LeftPar"){z <- z %>%
# mutate(`2014` = ifelse(`2013` != 1 & `2014` == 1, 1, 0),
# `2015` = ifelse(`2014` != 1 & `2015` == 1, 1, 0))
# }else{z <- z %>%
# mutate(`2015` = ifelse(`2013` != 1 & `2015` == 1, 1, 0))}
z <- z %>% mutate(
Event12 = ifelse(`2005` == 0 & rowSums(cbind(`2006`,`2007`,`2008`,`2009`),na.rm = T) > 0, 1, 0),
Event23 = ifelse(Event12 == 0 & rowSums(cbind(`2010`,`2011`,`2012`,`2013`),na.rm = T) > 0, 1, 0))
if(event != "LeftPar"){
z <- z %>% mutate(
Event3p = ifelse(rowSums(cbind(`2014`, `2015`), na.rm = T) > 0, 1, 0)) #%>%
# gather(key = le.group, value = le.value, `2006`:Event3p)
} else{
z <- z %>% mutate(Event3p = ifelse(`2015` >= 1, 1, 0)) #%>%
# gather(key = le.group, value = le.value, `2006`:Event3p)
}
}
# missing moving out of parental home and how often a child was born
le_dat <- dat %>% ungroup() %>%
separate(new_name, c("type", "Event"), sep = "_") %>%
filter(type == "LE" & #PROC_SID %in% unique(bfi_wide$PROC_SID) &
Event %in% c("Married", "Divorce", "MoveIn","SepPart", "PartDied",
"ChldMvOut", "ChldBrth", "MomDied", "DadDied", "ParDied",
"Unemploy", "Retire", "FrstJob", "LeftPar")) %>%
mutate(value = ifelse(Event == "Retire" | Event == "Unemploy", mapvalues(value, c(2,1), c(0,1)),
ifelse(Event == "FrstJob", mapvalues(value, seq(1,6,1), c(1,rep(0,5))),value))) %>%
group_by(Event) %>%
nest() %>%
mutate(event.dat = map2(data, Event, event_fun)) %>%
unnest(event.dat, .drop = T)
le_dat <- le_dat %>% select(Event:Dem_Sex, Event12, Event23) %>%
mutate(le.group = ifelse(Event12 == 1 | Event23 == 1, 1, 0),
le.group = ifelse(is.na(Event12) == T & is.na(Event23) == T, NA_real_, le.group))Clean Matching Data
The matching data is a little more complicated. Ultimately, we want to end up with 53 variables for matching, but these will be created from compositing more than one thousand responses collected between 1984 and 2005. To do this, we’ll again have to read in some raw data since some of the matching variables come from different datasets than the personality and life event data.
all.old.cols <- (meta %>% filter(class %in% c("proc") & Year == "0"))$Item
all.new.cols <- (meta %>% filter(class %in% c("proc") & Year == "0"))$new_name
# create short function to read in separate files for each wave
read_fun <- function(file, year){
print(year)
old.names <- (meta %>% filter(Year == year & class %in% c("match", "proc") & Include == "Yes"))$Item
new.names <- (meta %>% filter(Year == year & class %in% c("match", "proc") & Include == "Yes"))$new_name
z <- haven::read_sav(sprintf("%s/data/sav_files/%sp.sav", data_path, file)) %>%
left_join(haven::read_sav(sprintf("%s/data/sav_files/%spkal.sav", data_path, file))) %>%
left_join(haven::read_sav(sprintf("%s/data/sav_files/%sh.sav", data_path, file))) %>%
left_join(haven::read_sav(sprintf("%s/data/sav_files/%spequiv.sav", data_path, file))) %>%
left_join(haven::read_sav(sprintf("%s/data/sav_files/%spgen.sav", data_path, file))) %>%
left_join(haven::read_sav(sprintf("%s/data/sav_files/%shbrutto.sav", data_path, file))) %>%
# left_join(haven::read_sav(sprintf("%s/data/sav_files/%shost.sav", data_path, file))) %>%
# left_join(haven::read_sav(sprintf("%s/data/sav_files/%spost.sav", data_path, file))) %>%
select(one_of(all.old.cols), one_of(old.names)) %>%
setNames(c(all.new.cols, new.names)) %>%
mutate_all(funs(mapvalues(., seq(-1,-7,-1), c(NA, 0, rep(NA,5)), warn_missing = F))) %>%
group_by(PROC_SID) %>%
gather(key = new_name, value = value, -PROC_SID, -PROC_household, -Dem_DOB, -Dem_Sex) %>%
left_join(meta %>% filter(Year == year) %>% select(new_name, rev_code, mini, maxi, rule)) %>%
mutate(value = ifelse(rev_code == 1, reverse.code(keys=-1, items=value, mini=mini, maxi=maxi), value)) %>%
select(-rev_code, -mini, -maxi)
}
health.old.cols <- (meta %>% filter(dataset == "health"))$Item
health.new.cols <- (meta %>% filter(dataset == "health"))$new_name
health <- tbl_df(haven::read_sav(sprintf("%s/data/sav_files/health.sav", data_path))) %>%
filter(valid == 1 & syear < 2006) %>%
select(one_of(all.old.cols), one_of(health.old.cols)) %>%
setNames(c(all.new.cols, health.new.cols)) %>%
mutate_all(funs(mapvalues(., seq(-1,-7,-1), c(NA, 0, rep(NA,5)), warn_missing = F))) %>%
gather(key = new_name, value = value, -PROC_SID, -PROC_household, -Year) %>%
left_join(meta %>% filter(dataset == "health") %>% select(new_name, rule))
match.dat <- tibble(
Year = seq(1984, 2004,1),
file = c(letters[1:21])) %>%
mutate(data = map2(file, Year, read_fun)) %>%
unnest(data) %>%
full_join(health) %>%
mutate(Dem_DOB = ifelse(Dem_DOB < 1850, NA_real_, Dem_DOB),
Dem_DOB = recode(Dem_DOB, `0` = NA_real_),
Dem_Sex = recode(Dem_Sex, `0` = NA_real_)) %>%
group_by(PROC_SID) %>%
mutate(
Dem_DOB = max(Dem_DOB, na.rm = T),
Dem_Sex = max(Dem_Sex, na.rm = T)) #%>%Now, we need to recode the variables before compositing. In some cases, there are ordinal responses that make sense to collapse across levels. In other cases, we have to deal with missing data and such.
match.dat <- match.dat %>%
mutate(value = ifelse(new_name == "Psych_OthWorr" & value >= 1, 1, value),
value = ifelse(new_name %in% c("Bkgr_DadEdu", "Bkgr_MomEdu"),
mapvalues(value, c(0,6,7,1,2,9,3,4,5), rep(0:2, each = 3), warn_missing = F), value),
value = ifelse(new_name == "Bkgr_Edu" & value > 0, 1, value),
value = ifelse(new_name %in% c("Fnc_HouseAssist", "HH_Internet"),
recode(value, `2` = 0), value),
value = ifelse(new_name == "HH_CndHouse", mapvalues(value, c(2,3,4,1), c(0,0,0,1), warn_missing = F), value),
value = ifelse(new_name == "Bkgr_MarStat", mapvalues(value, c(2,6,7,1,3,4,5),
c(0,0,0,1,2,2,2), warn_missing = F), value),
value = ifelse(new_name == "HH_ClnHlp", mapvalues(value, c(3,1,2), c(0,1,1), warn_missing = F), value),
value = ifelse(new_name == "", mapvalues(value, c(98, 11,12,13,15,16, 21, 22, 31,32,33,34,
35,36,37,38, 14,41,44,42,43), c(0, rep(1, 15), rep(2,3), rep(3,2)), warn_missing = F), value),
Dem_DOB = ifelse(is.infinite(Dem_DOB) == T | is.nan(Dem_DOB) == T, NA, Dem_DOB),
Dem_Sex = ifelse(is.infinite(Dem_Sex) == T | is.nan(Dem_Sex) == T, NA, Dem_Sex))Now it’s time to composite some variables. To do so, we need to create a small function that will do this for us. We have 5 ways we’ll composite these variables: mean, mode, sum, selecting a specific year, and max. Which of these functions to use for different variables is indexed in our codebook, so all we have to do is use that as a guide for which to call.
After that’s done, we will join this matching data with teh BFI data for use in matching (for socialization).
# create a small function for calculating the mode
Mode <- function(x) {
ux <- unique(x)
ux <- ux[!is.na(ux)]
ux[which.max(tabulate(match(x, ux)))]
}
sum_fun <- function(df, Rule){
fun_call <- function(x, rule){
switch(rule,
average = mean(x, na.rm = T),
mode = Mode(x)[1],
sum = sum(x, na.rm = T),
select = unique(x)[1],
max = max(x, na.rm = T))
}
df %>%
group_by(PROC_SID, new_name, Dem_DOB, Dem_Sex, PROC_household) %>%
summarize(value = fun_call(value, Rule)) %>%
mutate(value = ifelse(is.nan(value) == T, NA,
ifelse(is.infinite(value) == T, NA, value))) %>%
ungroup()
}
match.dat.wide <- match.dat %>%
group_by(rule) %>%
nest() %>%
mutate(data = map2(data, rule, possibly(sum_fun, NA_real_))) %>%
unnest(data, .drop = T) %>%
select(-rule) %>% # get rid of the rule variables
spread(key = new_name, value = value) %>% # change to wide format
left_join(bfi_match)Match Subjects Across Data Categories
Our exclusion criteria is that participants had to have at least one wave of personality, matching, and life event data. So we’ll get a vector of the subjects for each type, then merge these into a vector of subjects who appear in all 3. Then we’ll trim out data sets for use in imputation, matching, and selection and socialization models for later.
bfi_subs <- unique(bfi_wide$PROC_SID)
match_subs <- unique(match.dat.wide$PROC_SID)
le_subs <- unique(le_dat$PROC_SID)
subs <- bfi_subs[bfi_subs %in% match_subs]
subs <- subs[subs %in% le_subs]
match.dat.wide <- match.dat.wide %>% filter(PROC_SID %in% subs)
bfi_wide <- bfi_wide %>% filter(PROC_SID %in% subs)
bfi_match <- bfi_match %>% filter(PROC_SID %in% subs)
le_dat <- le_dat %>% filter(PROC_SID %in% subs)
save(match.dat, match.dat.wide, bfi_wide, le_dat, bfi_match,
file = sprintf("%s/results/data.RData", data_path))Multiple Imputation
Missing Data Frame
First, we check the missingness patterns of the match data by converting it to a missing data frame class object using the missing_data.frame() function in the mi package in R. We then use the image() function to graphically depict the missingness patterns.
# MI doesn't like tibbles, so we need to unclass and reclass the data
match.dat.wide <- data.frame(unclass(match.dat.wide)) # mi doesn't like tibbles
mdf <- missing_data.frame(match.dat.wide)
pdf(sprintf("%s/plots/%s.pdf", data_path, "mdf"), width = 9, height = 6.5)
image(mdf)
dev.off()Now, we want to ensure that missing_data.frame() has correctly detected the type of variable (nominal, ordinal, etc.), so that missing data will be imputed using the correct link function.
des <- describe(match.dat.wide,fast = T)
mdf <- change(data = mdf, y = rownames(des)[des$range >= 3],
what = "type", to = "continuous")Multiple Imputation Procedure
Now, we use the mi() function in the mi package in to complete the multiple imputation procedure. We create 10 imputed data sets (by setting n.chains to 10) and use 20 iterations for each. We have a lot of variables and a lot of observations, so we use parallel processing to run the procedure.
mi.res <- mi(mdf, n.chains = 10, n.iter = 20, parallel = T)Now we compare the missingness patterns before and after imputation and see that we no longer have missing data.
pdf(sprintf("%s/plots/%s.pdf", data_path, "mi"), width = 9, height = 6.5)
image(mi.res)
dev.off()And grab the imputed data sets from the MI procedure using the complete() function from the mi package, which saves them in a list.
complete_fun <- function(mi){
clean_fun <- function(df){df %>% select(-contains("missing_"))}
tibble(chain = 1:10,
imp.data = mi::complete(mi)) %>%
mutate(imp.data = map(imp.data, clean_fun)) %>%
unnest(imp.data)
}
imp.data <- complete_fun(mi.res)Because we are estimating Bayesian growth curve models in STAN using the brms package, I’ll take the imputed data and grab the imputed personality variables. STAN isn’t cool with missing data, so this ensures that people won’t get thrown out.
bfi.imp <- unique(imp.data %>%
select(chain, PROC_SID, A1_1:O3_3) %>%
gather(key = item, value = value, A1_1:O3_3) %>%
separate(item, c("item", "wave"), sep = "_") %>%
separate(item, c("Trait", "item"),1) %>%
unite(item, wave, item, sep = "_") %>%
mutate(item = sprintf("T%s", item)) %>%
spread(key = item, value = value) %>%
left_join(bfi_wide %>% select(PROC_SID, sex12:age.c3)))Now, we’ll restructure the imputed data for the matching procedure. We’ll keep all the matching variables as well as imputed personality items from 2005. The personality items will only be used in matching for tests of socialization.
psw.imp.data <- imp.data %>%
gather(key = item, value = value, A1_1:O3_3) %>%
separate(item, c("item", "wave"), sep = "_") %>%
separate(item, c("Trait", "item"), 1) %>%
filter(wave == 1) %>%
group_by(chain, Trait, PROC_SID, wave) %>%
summarize(M = mean(value, na.rm = T)) %>%
ungroup() %>%
unite(Trait, Trait, wave, sep = "_") %>%
spread(key = Trait, value = M) %>%
full_join(imp.data %>% select(chain:Val_Chrch))
save(mdf, mi.res, file = paste(data_path, "results/mi_dat.RData", sep = "/"))
save(imp.data, bfi.imp, psw.imp.data, file = paste(data_path, "results/mi_dat_small.RData", sep = "/"))
rm("mi.res")We only need the columns that aren’t meant to define missingness patterns, so we write a simple function to do that for each element in the list of imputed data sets. Then we put them into a dataframe in which each of the imputed data sets is saved in a cell of the data frame, which will make it much easier to use for propensity score weighting and growth curve modeling.
load(url(paste(data_path, "results/mi_dat_small.RData", sep = "/")))
nested.psw <- crossing(
Event = unique(le_dat$Event),
match_set = c("selection", "socialization"),
chain = 1:10
)Propensity Score Matching
Then, we perform propensity score matching using our imputed datasets using the MatchIt package. We then extract the matched sets and merge in our predictor variables. To test the effectiveness of the propensity score matching procedure, we examine the average standardized effect size in the balance tables. Minimal effect sizes are candidates for beng dropped from the propensity score matching, and large effect sizes mean our matching procedure wasn’t effective. We can also examine these using balance plots.
The function below looks a little complicated. But here’s what it does. We have to run the matching procedure for each of the 10 chains of imputed data for each life event for both socialization (matching variables + personality) and selection (matching variables only) sets. Some events are harder to match than others because the sample sizes are more unbalanced. For these, we’ll use a larger matching ratio (8:1 instead of our default 4:1). In addition, some sets remained unbalanced despite teh change in the ratio. For these sets, we also reduced the caliper width from .25\(\sigma\) to .05\(\sigma\) (Partner Separation, Unemployment, Retirement, and First Job) or .01\(\sigma\)(Move in with a Partner). We do this because it is critical for our tests of socialization and selection that our samples be matched.
# this function actually runs the propensity score weighting procedure
psw_fun <- function(event, Chain, match_set){
print(sprintf("%s Chain %s", event, Chain))
Ratio <- ifelse(event %in% c("ChldMvOut", "ParDied", "MoveIn"), 8, 4)
Caliper <- ifelse(match_set == "socialization" &
event %in% c("SepPart", "Unemploy", "Retire", "FrstJob"), .05,
ifelse(match_set == "socialization" &
event %in% c("MoveIn"), .01, .25))
df <- psw.imp.data %>% filter(chain == Chain) %>%
full_join(le_dat %>% filter(Event == event) %>%
select(Event, PROC_SID, le.group)) %>%
select(-chain, -Event)
if(event == "Retire"){df <- df %>% filter(2005 - Dem_DOB > 40)}
if(event == "FrstJob"){df <- df %>% filter(2005 - Dem_DOB < 40)}
df <- df[complete.cases(df),]
df <- data.frame(unclass(df))
if(match_set == "socialization"){to.match <- colnames(df)[-which(colnames(df) %in% c("PROC_SID","le.group",
paste(rep(c("A", "C", "E", "N", "O"), each = 2), rep(2:3, times = 5), sep = "_")))]}
else {to.match <- colnames(df)[-which(colnames(df) %in% c("PROC_SID","le.group",
paste(rep(c("A", "C", "E", "N", "O"), each = 3), rep(1:3, times = 5), sep = "_")))]}
match.formula <- as.formula(paste("le.group ~ ", paste(to.match, collapse=" + "), sep = " "))
y <- matchit(match.formula, data = df, method = "nearest", ratio = Ratio, caliper = Caliper)
}
# changing the data fed into psw into a data frame because it won't work with tibbles
psw_df <- function(psw){psw$data <- data.frame(psw$data); psw}
psw_df <- function(psw){
data.frame(match.data(psw))
}This function creates the balance table of the psw weights and filters the results into variables the matching procedure did not fix and those that it did.
unbalanced_fun <- function(x){
y <- summary(x, standardize = T)
raw <- y$sum.all %>%
mutate(var = rownames(.)) %>%
select(var, `Means Treated`, `Means Control`, `Std. Mean Diff.`)
smalldiff.var <- raw %>% filter(abs(`Std. Mean Diff.`) <= .05)
matched <- y$sum.matched %>%
mutate(var = rownames(.)) %>%
select(var, `Means Treated`, `Means Control`, `Std. Mean Diff.`)
unbalanced.var <- matched %>% filter(abs(`Std. Mean Diff.`) >= .2)
return(list(raw = raw, matched = matched,
unbalanced = unbalanced.var,smalldiff = smalldiff.var))
}Let’s run it now!
nested.psw <- nested.psw %>%
filter((chain == 1 & match_set == "socialization") |
match_set == "selection") %>%
mutate(psw = pmap(list(Event, chain, match_set), psw_fun),
psw.df = map(psw, possibly(psw_df, NA_real_)),
bal.df = map(psw, unbalanced_fun),
raw = map(bal.df, ~.$raw),
matched = map(bal.df, ~.$matched),
unbal.tab = map(bal.df, possibly(~.$unbalanced, NA_real_)),
smalldiff.tab = map(bal.df, possibly(~.$smalldiff, NA_real_)))
save(nested.psw, file = paste(data_path, "results/psw.RData", sep = "/"))
nested.psw <- nested.psw %>% select(-psw)
save(nested.psw, file = paste(data_path, "results/psw_small.RData", sep = "/"))Balance Plots
In these plots, substantial reductions in effect sizes are observed for most variables (blue lines), with only one variable showing an increase in effect size (red lines), but only a seemingly trivial increase. Closed red circles indicate a statistically significant difference, many of which occur before weighting, none after.
plot_fun <- function(df, event, set){
plot <- df %>%
mutate(type = factor(type, level = c("raw", "matched"))) %>%
ggplot(aes(x = type, y = `Std. Mean Diff.`)) +
scale_y_continuous(limits = c(-5,5), breaks = seq(-5,5,2)) +
geom_point() +
geom_line(aes(group = var), size = .25, alpha = .8) +
labs(x = NULL, y = "Standardized Mean Difference",
title = sprintf("%s", event)) +
facet_wrap(~chain, ncol = 2) +
theme_classic()
# ggsave(sprintf("%s/plots/psw_bal_%s.png", data_path, event), width = 8, height = 10)
}
print_plot_fun <- function(p, Chain){
p$main <- sprintf("Imputed dataset %s", gsub("chain.", "", Chain))
p
}
nested.psw %>%
unnest(raw) %>% mutate(type = "raw") %>%
full_join(nested.psw %>% unnest(matched) %>% mutate(type = "matched"))%>%
group_by(Event) %>%
nest() %>%
mutate(plot = map2(data, Event, plot_fun))
par(mfrow = c(2,5))
nested.psw <- nested.psw %>%
mutate(plots = map(psw, possibly(~plot(., plots = "es"), NA_real_)),
plots = map2(plots, chain, possibly(print_plot_fun, NA_real_))) Balance Tables
Once propensity scores are estimated, bal.table() produces a table that shows how well the resulting weights succeed in manipulating the groups so that they match on pre-adolescent matching characteristics.
After matching, none of the variables remained unbalanced (standardized mean differences > .2), so there is no table of unbalanced variables (yay!).
load(url(paste(data_path, "results/psw_small.RData", sep = "/")))
load(url(paste(data_path, "results/mi_dat_small.RData", sep = "/")))
# this table shows variables that are not matched after weighting
# unbal.tab <- nested.psw %>%
# filter(match_set == "socialization") %>%
# unnest(unbal.tab, .drop = T) %>%
# group_by(Event, var) %>%
# #summarize_at(vars(`Means Treated`:`Std. Mean Diff.`), funs(mean(., na.rm = T)))
# summarize(mean = mean(`Std. Mean Diff.`, na.rm = T)) %>%
# spread(key = Event, value = mean)
# kable(unbal.tab, "html", longtable = T, booktabs = T, digits = 2,
# caption = "Unbalanced Variables after Propensity Score Weighting") %>%
# kable_styling(bootstrap_options = c("striped","repeat_header"),full_width = F)
# this table shows variables that were already matched prior to weighting
smalldiff.tab <- nested.psw %>%
unnest(smalldiff.tab, .drop = T) %>%
group_by(Event, var) %>%
#summarize_at(vars(`Means Treated`:`Std. Mean Diff.`), funs(mean(., na.rm = T)))
summarize(mean = mean(`Std. Mean Diff.`, na.rm = T)) %>%
spread(key = Event, value = mean)
kable(smalldiff.tab, "html", longtable = T, booktabs = T, digits = 2,
caption = "Balanced Variables after Propensity Score Weighting") %>%
kable_styling(bootstrap_options = c("striped","repeat_header"),full_width = F)| var | ChldBrth | ChldMvOut | DadDied | Divorce | FrstJob | LeftPar | Married | MomDied | MoveIn | ParDied | PartDied | Retire | SepPart | Unemploy |
|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|
| A_1 | 0.01 | 0.01 | 0.02 | -0.04 | -0.04 | |||||||||
| Act_Volunteer | 0.04 | -0.04 | 0.02 | 0.02 | ||||||||||
| Bkgr_DadPres1 | -0.02 | 0.02 | ||||||||||||
| Bkgr_DisabStat1 | -0.05 | -0.03 | ||||||||||||
| Bkgr_Edu1 | 0.04 | -0.02 | ||||||||||||
| Bkgr_MarStat.L | -0.04 | |||||||||||||
| Bkgr_MarStat.Q | 0.00 | |||||||||||||
| Bkgr_PGovIncome | 0.00 | 0.04 | ||||||||||||
| Bkgr_UrbOrRur2 | -0.03 | -0.01 | 0.02 | 0.00 | -0.01 | 0.04 | -0.03 | 0.00 | -0.02 | |||||
| C_1 | 0.00 | 0.00 | 0.04 | -0.04 | ||||||||||
| Dem_DOB | 0.00 | |||||||||||||
| Dem_Sex1 | -0.02 | -0.01 | -0.01 | 0.04 | 0.03 | 0.04 | -0.02 | 0.01 | 0.01 | 0.02 | ||||
| Dem_Sex2 | 0.02 | 0.01 | 0.01 | -0.04 | -0.03 | -0.04 | 0.02 | -0.01 | -0.01 | -0.02 | ||||
| E_1 | 0.05 | -0.01 | -0.01 | 0.00 | 0.00 | 0.04 | ||||||||
| Fnc_HouseAssist1 | 0.00 | -0.03 | ||||||||||||
| Fnc_StudGrnt1 | -0.03 | 0.02 | ||||||||||||
| Fnc_UnempBen | -0.04 | 0.01 | 0.03 | -0.01 | 0.01 | -0.01 | ||||||||
| HH_BrothPres1 | -0.03 | |||||||||||||
| HH_ClnHlp1 | 0.02 | 0.01 | -0.05 | 0.04 | 0.01 | 0.02 | ||||||||
| HH_CndHouse1 | 0.02 | -0.04 | ||||||||||||
| HH_ColTV1 | 0.01 | 0.04 | -0.04 | -0.04 | -0.02 | 0.01 | 0.01 | |||||||
| HH_NumPer | 0.01 | |||||||||||||
| HH_NumPer15to18 | 0.03 | 0.05 | 0.00 | 0.03 | 0.04 | |||||||||
| HH_NumPerBel14 | -0.04 | |||||||||||||
| HH_SisPres1 | 0.00 | |||||||||||||
| Hlth_BMI | 0.01 | |||||||||||||
| Hlth_BodPain | 0.02 | -0.05 | 0.00 | |||||||||||
| Hlth_EmoRole | -0.01 | -0.04 | 0.02 | -0.05 | 0.04 | -0.05 | -0.03 | -0.02 | ||||||
| Hlth_GenHlth | 0.01 | -0.03 | ||||||||||||
| Hlth_HeightCM | -0.03 | 0.02 | -0.02 | 0.01 | ||||||||||
| Hlth_HlthInsr | -0.03 | 0.01 | -0.04 | -0.03 | ||||||||||
| Hlth_MntlHlth | -0.05 | |||||||||||||
| Hlth_NumDrVisits | 0.04 | -0.03 | ||||||||||||
| Hlth_PhysFunc | -0.04 | |||||||||||||
| Hlth_PhysHlth | -0.02 | |||||||||||||
| Hlth_PhysProb | 0.01 | -0.02 | ||||||||||||
| Hlth_PhysRole | 0.03 | 0.00 | ||||||||||||
| Hlth_SocFunc | 0.00 | 0.04 | 0.00 | -0.02 | -0.02 | |||||||||
| Hlth_Vitality | -0.01 | 0.03 | -0.01 | -0.02 | 0.01 | 0.03 | ||||||||
| Hlth_WeightKG | -0.02 | 0.05 | 0.04 | 0.04 | ||||||||||
| N_1 | -0.01 | -0.02 | 0.01 | 0.04 | -0.04 | 0.03 | -0.03 | |||||||
| O_1 | -0.01 | 0.03 | 0.03 | 0.00 | 0.02 | -0.03 | ||||||||
| PROC_household | 0.03 | -0.05 | 0.00 | 0.05 | 0.02 | -0.01 | ||||||||
| Psych_LifeSat | -0.05 | -0.02 | -0.01 | -0.01 | -0.04 | 0.04 | ||||||||
| Psych_OthWorr1 | 0.05 | 0.02 | -0.05 | 0.02 | 0.02 | 0.04 | -0.02 | -0.01 | 0.00 | |||||
| Psych_SatFam | 0.04 | 0.04 | 0.04 | -0.02 | 0.02 | 0.00 | -0.01 | |||||||
| Psych_SatHealth | 0.01 | -0.04 | 0.03 | |||||||||||
| Psych_SatIncome | -0.03 | 0.04 | 0.01 | |||||||||||
| Psych_WorrCrm | 0.03 | -0.04 | -0.05 | 0.01 | ||||||||||
| Rel_RelDad | 0.00 | 0.02 | ||||||||||||
| Rel_RelMom | 0.04 | |||||||||||||
| Soc_SocGath | 0.04 | 0.02 | ||||||||||||
| Soc_VisFam | 0.00 | 0.03 | 0.05 | 0.02 | -0.03 | |||||||||
| Soc_VisNghbr | 0.03 | -0.04 | ||||||||||||
| Val_Chrch | -0.02 | 0.03 | -0.03 |
Plots
In addition to the balance plots that highlight the attenuation of mean differences in general, it is also interesting to see this for specific variables. To do that, we’ll calculate Cohen’s d (standardized mean differences) between groups for each matching variable in both the matched and unmatched sets. Then, we’ll plot them together.
Standardized Mean Differences
cohens_d <- function(x, y) {
lx <- length(x)- 1; ly <- length(y)- 1
md <- mean(x, na.rm = T) - mean(y, na.rm = T) ## mean difference (numerator)
csd <- lx * var(x, na.rm = T) + ly * var(y, na.rm = T)
csd <- csd/(lx + ly); csd <- sqrt(csd) ## common sd computation
cd <- md/csd ## cohen's d
return(cd)
}
d_fun <- function(df, Var, chain){
dat <- df %>% filter(chain == chain)
d <- with(dat, cohens_d(value[le.group == 1], value[le.group == 0]))
}
levs <- c(colnames(psw.imp.data)[-c(1,2,10)], "Age")
diff <- psw.imp.data %>% filter(chain == 1) %>%
mutate(match_set = "Unmatched") %>%
full_join(le_dat %>% select(PROC_SID, Event, le.group)) %>%
# filter(Event %in% c("FrstJob", "MoveIn", "PartDied"))) %>%
filter(!is.na(Event) & !is.na(match_set)) %>%
full_join(nested.psw %>%
filter(match_set == "socialization" & chain == 1) %>%
# filter(Event %in% c("FrstJob", "MoveIn", "PartDied")) %>%
unnest(psw.df) %>%
select(-weights)) %>%
mutate(Dem_Sex = as.numeric(as.character(Dem_Sex)),
Age = 2005-Dem_DOB) %>% select(-Dem_DOB) %>%
rename(BFI_E.W1 = E_1, BFI_A.W1 = A_1, BFI_C.W1 = C_1, BFI_N.W1 = N_1, BFI_O.W1 = O_1) %>%
mutate_if(is.factor, funs(as.numeric(as.character(.)))) %>%
gather(key = var, value = value, -chain, -PROC_SID, -Event, -PROC_household, -match_set, -le.group) %>%
separate(var, c("Category", "Item"), sep= "_", remove = F) %>%
# mutate(var = factor(var, levels = rev(levs))) %>%
group_by(Event, match_set, var, chain, Category, Item) %>%
nest() %>%
mutate(d = pmap(list(data, var, chain), d_fun))
# save(diff, file = sprintf("%s/results/matching_diag.RData", data_path))
# save(diff, file = sprintf("%s/results/mean_diff.RData", data_path))
events <- tibble(
old = c("none", "Married", "MoveIn", "Divorce", "SepPart", "PartDied", "LeftPar",
"ChldMvOut", "ChldBrth", "ParDied", "Unemploy", "Retire", "FrstJob"),
new = c("Mean", "Marriage", "Moved in with Partner", "Divorce", "Separation from Partner",
"Death of Partner/Spouse", "Leaving Parental Home", "Child Leaves Home",
"Birth of Child", "Death of Parent", "Unemployment", "Retirement", "First Job"),
breaks = c("Mean", "Marriage", "Moved in\nwith Partner", "Divorce", "Separation\nfrom Partner",
"Death of\nPartner/Spouse", "Leaving\nParental Home", "Child Leaves\nHome",
"Birth of\nChild", "Death of\nParent", "Unemployment", "Retirement", "First Job")
)
diff_fun <- function(event){
p<-diff %>% unnest(d, .drop = T) %>% #filter(Category == "HH") %>%
filter(var != "distance") %>%
filter(Event == event) %>%
group_by(Event, match_set, var) %>%
summarize(d = mean(d)) %>% ungroup() %>%
mutate(match_set = recode(match_set, `socialization` = "Matched"),
Event = mapvalues(Event, events$old, events$breaks),
Event = factor(Event, levels = events$breaks)) %>%
ggplot(aes(x = var, y = d, shape = match_set)) +
scale_shape_manual(values = c(19,1)) +
scale_y_continuous(limits = c(-1.5, 1.5), breaks = seq(-1, 1, 1)) +
geom_hline(aes(yintercept = 0), linetype = "dashed", size = .25) +
geom_point(size = 1.5) +
labs(y = "Cohen's d", x = NULL, shape = NULL, title = event) +
# coord_flip() +
facet_grid(.~Event) +
theme_classic() +
theme(legend.position = "bottom",
axis.text.x = element_text(face = "bold", size = rel(.7), angle = 45, hjust = 1),
axis.text.y = element_text(face = "bold", size = rel(1.2)),
axis.title = element_text(face = "bold", size = rel(1.2)),
strip.text = element_text(face = "bold"),
plot.title = element_text(face = "bold", size = rel(1.2), hjust = .5),
legend.text = element_text(face = "bold"),
legend.title = element_text(face = "bold", size = rel(1.2)))
print(p)
}for (x in events$old[events$old != "none"]){
cat('\n\n####', x, '\n\n ')
diff_fun(x)
}Married
MoveIn
Divorce
SepPart
PartDied
LeftPar
ChldMvOut
ChldBrth
ParDied
Unemploy
Retire
FrstJob
# ggsave(sprintf("%s/results/plots/psm_match.png", data_path), width = 10, height = 4)Propensity Scores
What about the distribution of propensity scores? Well, we can also look at the propensity score matched sets. The chunk below will take the propensity scores (“distances” from the matchit object) and plot those as a function of group and weight. Larger points on the graph indicate participants who were weighted more. For example, a parpticipant with a from the Event group who had four matches from the no event group would be weighted 1, while each of their matches would be weighted .25. We don’t actually use the weights because are matching, but they are useful to visualize essentially how the matching procedure approximates using survey weights.
sample_fun <- function(df){
df %>% filter(row_number() %in% sample(1:nrow(df), 200))
}
nested.psw %>%
filter(!Event %in% c("MomDied", "DadDied")) %>%
mutate(sample = map(psw.df, sample_fun)) %>%
unnest(sample) %>%
select(Event:PROC_SID, le.group, distance, weights) %>%
mutate(le.group = mapvalues(le.group, c(0,1), c("No Event", "Event")),
Event = mapvalues(Event, events$old, events$breaks),
Event = factor(Event, levels = events$breaks)) %>%
filter(chain == 1 & match_set == "socialization") %>%
ggplot(aes(x = le.group, y = distance)) +
scale_size_continuous(range = c(.1,3)) +
scale_y_continuous(breaks = seq(0,1,.5), labels = c("0", ".5", "1")) +
geom_jitter(aes(size = distance), shape = 1, width = .15) +
labs(x = NULL, y = "Propensity Score") +
coord_flip() +
facet_wrap(~Event, nrow = 3) +
theme_classic() +
theme(legend.position = "none",
strip.text = element_text(face = "bold", size = rel(1)),
axis.text = element_text(face = "bold", size = rel(1.2)),
axis.title = element_text(face = "bold", size = rel(1.2)))
# ggsave(file = sprintf("%s/results/plots/ps_plot.png", data_path),
# width = 10, height = 6)Group Differences
Finally, let’s make a table of the some of the demographic and sample characteristics of the matched samples. Below, I create a table that includes the sample sizes of each group for each event for both matched (socialization) and raw samples, mean and standard deviation of age of people who experienced the events and the percent of people who experienced the events who were women.
events <- tibble(
old = c("none", "Married", "MoveIn", "Divorce", "SepPart", "PartDied", "LeftPar",
"ChldMvOut", "ChldBrth", "ParDied", "Unemploy", "Retire", "FrstJob"),
new = c("Mean", "Marriage", "Moved in with Partner", "Divorce", "Separation from Partner",
"Death of Partner/Spouse", "Leaving Parental Home", "Child Leaves Home",
"Birth of Child", "Death of Parent", "Unemployment", "Retirement", "First Job"),
breaks = c("Mean", "Marriage", "Moved in\nwith Partner", "Divorce", "Separation\nfrom Partner",
"Death of\nPartner/Spouse", "Leaving\nParental Home", "Child Leaves\nHome",
"Birth of\nChild", "Death of\nParent", "Unemployment", "Retirement", "First Job")
)
options(knitr.kable.NA = '')
options(papaja.na_string = "")
options(knitr.table.format = "html")
dem <- le_dat %>% mutate(chain = 1, data = "raw", Dem_Sex = factor(Dem_Sex)) %>%
full_join(nested.psw %>%
filter(match_set == "socialization" & chain == 1) %>%
unnest(psw.df) %>%
select(Event, chain, PROC_SID, Dem_DOB, Dem_Sex, le.group) %>%
mutate(data = "matched")) %>%
mutate(age = 2005-Dem_DOB) %>%
#select(-le.value)) %>%
filter(!is.na(le.group)) %>%
group_by(Event, data, chain) %>%
mutate(m.age = mean(age[le.group == 1], na.rm = T),
sd.age = sd(age[le.group == 1], na.rm = T))
dem %>%
group_by(Event, chain, data, le.group, m.age, sd.age) %>%
dplyr::summarize(n = n()) %>%
full_join(dem %>% group_by(Event, data, Dem_Sex, le.group, chain) %>% summarize(n = n()) %>%
unite(temp, Dem_Sex, le.group, sep = ".") %>% spread(key = temp, value = n) %>%
mutate(perc_women = `2.1`/rowSums(cbind(`1.1`, `2.1`))*100) %>%
select(Event, chain, data, perc_women)) %>%
group_by(Event, data, le.group) %>%
summarize_at(vars(m.age:perc_women), funs(mean(., na.rm = T))) %>%
ungroup() %>%
spread(key = le.group, value = n) %>%
mutate(Frequency = sprintf("%.0f (%.0f)", `1`, (`1` + `0`))) %>%
mutate_at(vars(m.age:perc_women), funs(sprintf("%.2f", .))) %>%
select(-`0`, -`1`) %>%
gather(key = measure, value = value, m.age:Frequency) %>%
filter(data == "matched") %>%
unite(data, measure, data, sep = ".") %>%
spread(key = data, value = value) %>%
select(Event, contains("Frequency"), contains("age"), everything()) %>%
filter(!(Event %in% c("MomDied", "DadDied"))) %>%
mutate(Event = mapvalues(Event, events$old, events$new),
Event = factor(Event, levels = events$new)) %>%
arrange(Event) %>%
kable(., "html", booktabs = T, escape = F,
col.names = c("Life Event", rep(c("Matched"), times = 4)),
caption = "Descriptive Statistics of Individuals Who Experienced a Specific Major Life Event") %>%
kable_styling(full_width = F, font_size = 8) %>%
#kable_styling(full_width = F) %>%
add_header_above(c(" " = 1, "Frequency" = 1, "M" = 1, "SD" = 1, "% women" = 1)) %>%
add_header_above(c(" " = 2, "Age in 2005" = 2, " " = 1)) %>%
footnote(general = "M = mean age in 2005; SD = standard deviation of age in 2005.") | Life Event | Matched | Matched | Matched | Matched |
|---|---|---|---|---|
| Marriage | 965 (4458) | 33.36 | 11.04 | 51.40 |
| Moved in with Partner | 884 (4890) | 31.56 | 10.77 | 52.94 |
| Divorce | 362 (1803) | 40.73 | 8.43 | 58.84 |
| Separation from Partner | 947 (4531) | 35.62 | 11.36 | 56.60 |
| Death of Partner/Spouse | 390 (1927) | 64.08 | 11.74 | 70.26 |
| Leaving Parental Home | 374 (1528) | 23.88 | 7.88 | 50.80 |
| Child Leaves Home | 1814 (7779) | 48.37 | 8.11 | 55.35 |
| Birth of Child | 1116 (4351) | 31.04 | 7.73 | 53.85 |
| Death of Parent | 1633 (10630) | 46.78 | 10.63 | 52.54 |
| Unemployment | 1910 (7506) | 38.45 | 12.68 | 51.94 |
| Retirement | 1166 (3221) | 60.55 | 10.79 | 54.03 |
| First Job | 374 (1122) | 22.11 | 3.85 | 52.67 |
| Note: | ||||
| M = mean age in 2005; SD = standard deviation of age in 2005. |
Selection Effects: Logistic Regressions
Made it to selection! For the selection models, we are using Bayesian logistic regression models to estimate the log odds of a participant experiencing an event as predicted by their personality in 2005. We’ll do this for each combination of trait, event, and chain (from the imputed data sets) with the matched sets as well as for each trait and event combination for the raw data from 2005. Doing both unmatched and matched will allow us to examine whether matching attenuates selection effects.
## Run Models
load(url(sprintf("%s/results/psw_small.RData", data_path)))
load(url(sprintf("%s/results/mi_dat_small.RData", data_path)))
bfi_selection <- crossing(
Event = c(unique(nested.psw$Event)),
Trait = unique(bfi_wide$Trait),
match = c("Matched", "Unmatched")
)
sel_fun <- function (event, trait, match){
k <- 1#ifelse(event == "ParDied" & trait == "A" & match == "Matched", 4, 1)
lapply(k:10, function(x){
print(paste(event, trait, match, x), sep = " ")
if(match == "Matched"){
subs <- unique((nested.psw %>%
filter(match_set == "selection" & chain == 1 & Event == event) %>%
unnest(psw.df))$PROC_SID)
df <- bfi.imp %>% filter(chain == x & Trait == trait & PROC_SID %in% subs)
} else {
df <- bfi.imp %>% filter(chain == x & Trait == trait)
}
df <- df %>%
mutate(value = rowMeans(cbind(T1_1, T1_2, T1_3), na.rm = T)) %>%
select(PROC_SID, sex12:age.c3, value) %>%
left_join(le_dat %>% filter(Event == event) %>% select(PROC_SID, le.group))
mod <- rstanarm::stan_glm(
le.group ~ value, family = binomial(link = "logit"), data = df
)
file <- sprintf("%s/results/selection/%s_%s_%s_chain%s.RData", data_path, trait, event, match, x)
save(mod, file = file)
rm(mod)
return(T)
})
}
start <- Sys.time()
bfi_selection <- bfi_selection %>%
mutate(b.sel.mod = pmap(list(Event, Trait, match), sel_fun))
end <- Sys.time()Pool Results
Because we had 10 imputations of the data, we have to pool the models across chains. Because these logistic regressions are not unwieldly in size, we can simply combine the samples from the Bayesian models across imputations and then save the result as a stan model object. Then, we can use default functions from the rstanarm package to extract the results.
stantab_fun <- function(event, trait, match){
models <- lapply(1:10, function(x){
print(x)
file <- sprintf("%s/results/selection/%s_%s_%s_chain%s.RData",
data_path, trait, event, match, x)
load(url(file))
return(mod)
})
sflist <- lapply(models, "[[", "stanfit")
models[[1]]$stanfit <- rstan::sflist2stanfit(sflist)
models[[1]]
}Table Results
Now we’ll run the pooling and extract teh samples (to visualize the posterior distributions).
bfi_selection <- crossing(
Event = c(unique(le_dat$Event)),
Trait = unique(bfi_wide$Trait),
match = c("Matched", "Unmatched")
) %>%
mutate(pool = pmap(list(Event, Trait, match), stantab_fun),
samples = map(pool, ~gather_samples(., value)),
qi = map(samples, mean_qi))
save(bfi_selection, file = sprintf("%s/results/selection.RData", data_path))Plots
Now, the fun part. Next, we’ll use the samples to plot the posterior distributions (after converting log odds to odds) as well as the mean estimates. We’ll do this for all trait x event x matching combinations. Note that within a Bayesian framework, we aren’t talking about significance, but I will highlight those models whose mean estimates suggest that in 95% of our samples, the highest posterior density interval does not include 1 (no difference in odds).
load(url(sprintf("%s/results/selection_small.RData", data_path)))
bfi_samples <- bfi_selection %>% unnest(samples) %>%
filter(term == "value") %>%
mutate(estimate = exp(estimate))
save(bfi_samples, files = sprintf("%s/results/selection_samples.RData", data_path))
bfi_selection <- bfi_selection %>% select(-samples, -pool)
save(sel.tab, bfi_selection, file = sprintf("%s/results/sel.tab.RData", data_path))load(url(sprintf("%s/results/selection_small.RData", data_path)))
load(url(sprintf("%s/results/selection_samples.RData", data_path)))
load(url(sprintf("%s/results/sel.tab.RData", data_path)))
bfi_qi <- bfi_selection %>% unnest(qi) %>%
filter(term == "value") %>%
mutate(sig = ifelse(sign(conf.low) == sign(conf.high), "sig", "ns")) %>%
mutate_at(vars(estimate:conf.high), funs(exp))
sel_plot_fun <- function(match_set){
p <- bfi_samples %>% full_join(bfi_qi %>% select(Event:term, sig)) %>%
filter(!(Event %in% c("MomDied", "DadDied"))) %>%
mutate(Trait = factor(Trait, levels = c("E", "A", "C", "N", "O")),
Event = mapvalues(Event, events$old, events$new),
Event = factor(Event, rev(events$new))) %>%
ggplot(aes(y = Event, x = estimate)) +
scale_x_continuous(limits = c(.4, 1.6), breaks = seq(.5, 1.5, .5), labels = c(".5", "1", "1.5")) +
# geom_density_ridges(data = . %>% filter(match == "Matched"), aes(fill = sig), rel_min_height = 0.025, scale = 4) +
geom_halfeyeh(data = . %>% filter(match == match_set), aes(fill = sig)) +
geom_blank(data = bfi_samples %>% filter(!(Event %in% c("MomDied", "DadDied"))) %>%
mutate(Trait = factor(Trait, levels = c("E", "A", "C", "N", "O")),
Event = mapvalues(Event, events$old, events$new),
Event = factor(Event, rev(events$new)))) +
# geom_halfeyeh(data = . %>% filter(match == "Matched" & sig == "sig"), fill = "lightgoldenrod") +
geom_vline(aes(xintercept = 1), linetype = "dashed") +
scale_fill_manual(values = c("gray", "lightgoldenrod")) +
labs(x = "OR", y = NULL, fill = NULL, color = NULL, title = match_set) +
facet_wrap(~Trait, nrow = 1) +
theme_classic() +
theme(legend.position = "none",
axis.text = element_text(face = "bold", size = rel(1.2)),
strip.text = element_text(face = "bold", size = rel(2)),
axis.title = element_text(face = "bold", size = rel(1.2)),
legend.text = element_text(face = "bold", size = rel(1.2)))
# ggsave(p, file = sprintf("%s/results/plots/selection_%s_posterior.png", data_path, match_set),
# width = 14, height = 6)
print(p)
}Sample Plots
Here is a small sample plot of the posterior distributions for use in my presnetation.
bfi_samples %>% full_join(bfi_qi %>% select(Event:term, sig)) %>%
filter(Trait %in% c("C", "E") & Event == "FrstJob") %>%
ggplot(aes(y = Event, x = estimate)) +
scale_x_continuous(limits = c(.4, 1.6), breaks = seq(.5, 1.5, .5),
labels = c(".5", "1", "1.5")) +
geom_halfeyeh(data = . %>% filter(match == "Matched"), aes(fill = sig)) +
geom_vline(aes(xintercept = 1), linetype = "dashed") +
scale_fill_manual(values = c("gray", "lightgoldenrod")) +
labs(x = "OR", y = NULL, fill = NULL, color = NULL) +
facet_wrap(~Trait, nrow = 1) +
theme_classic() +
theme(legend.position = "none",
axis.text = element_text(face = "bold", size = rel(1.2)),
strip.text = element_text(face = "bold", size = rel(2)),
axis.title = element_text(face = "bold", size = rel(1.2)),
legend.text = element_text(face = "bold", size = rel(1.2)))
# ggsave(file = sprintf("%s/results/plots/selection_example_plot.png", data_path),
# width = 6, height = 2)for (x in c("Matched", "Unmatched")){
cat(' \n\n### ', x, ' Selection Effects\n\n ', sep ="")
sel_plot_fun(x)
}Matched Selection Effects
Unmatched Selection Effects
Socialization Effects: Growth Models
The growth models are a little more complicated. Each of the 700 models themselves are so large after we run them (each is several GB), we can’t pool them like we did the Bayesian logistic regression models. We’d need a computer with some 30 GB of RAM and several additional months. So instead, we’re going to use Rubin’s Rules to pool the estimates. To do that, though, we need to extract certain elements from the models (e.g. standard errors of simple effects estimates, hypothesis tests of means, predicted values of both fixed and random effects), which is what the functions below do.
Define Functions
Standard Errors
Below we extract the variance-covariance matrix from each model and use it to estimate the standard error of the simple effects for the fixed effect parameters.
Hypothesis Tests
We’ve coded our growth curve models such that the coefficients tell us about mean differences among groups, but we are also interested in the mean estimates for each group because it’s useful to be able to say, for example, that (1) the groups changed differently such that (2) people who experienced the event didn’t cahnge (arrested development) while (3) people who didn’t experience the event increased in a trait over time. We’ll use the
hypothesis()function frombrms.Predicted Values
For plotting fixed and random trajectories over time, we want to get the predicted values from the model across the study period. To do this, we’ll use the
fitted()function for the fixed effects withre_formulaset to NA to ignore the random effects (for now). Then, we’ll create a data frame for each wave x event x wave x participant combination and use thepredict()function to get person level predictions.Pooling
Sorry in advance because there is a lot happening with this pooling function. The basic logic of it is this. These growth models are so large that most computers can only handle one at a time. So I need to load each model into R separately, extract everything I could possibly need from it, discard the model, and then save the stuff I need. I do this for each trait x event combination 10 times then pool my extracted results across those imputations using Rubin’s rules.
brm Models
The function below actually runs the
brmsmodels using thebrm()function. Note that the model formula is different for the life event v. no life event models. Because data storage is a problem. For each iteration of the models, we’ll filter out the responses we need. I don’t want to store these in the already unwieldily large model data frame. Note that we use the defaultbrmspriors in most cases, with the exception of the variance terms. For those we use a half-cauchy distribution, which is better than a student’s t because it doesn’t include negative terms (and variances can’t be negative.). In addition, some models had convergence issues. For these models, we’ll up the number of iterations from 2000 to 8000, the warmup from 1000 iterations to 4000 iterations, theadapt_deltato .99 rather than .8, and the maximum treedepth from 10 to 20. We’ll save each model to a .RData file that we’ll load in and pool using our other functions.Run Models
Finaly, let’s run the models. Please don’t actually do this unless you have months of your life or a super computer.
Pool Results
Now time to pool the results. Note that this will call the
pool_fun()I defined above that will load in and pool all the models. This will then save the results to yet another .RData file, that we’ll load in after pooling, which takes at least a week.Load Results
The models take days to weeks to run, and pooling takes at least an additional day, so we’re going to load in the results of the models, rather than the models themselves.
Table Results
Let’s create some tables. Figures are more fun, but sometimes it’s nice to see different values.
Intercepts
The table below shows the intercepts of those who did or didn’t experience events, along with the highest posterior density interval of that intercept and the standardized mean difference between groups in 2005.
Group Differences in Intercepts
The table below shows the differnces in intercepts between the group. Note that the matching procedure should have (and did) eliminate almost all of these.
[-0.01, 0.14]
[-0.09, 0.04]
[-0.07, 0.06]
[-0.13, 0.04]
[-0.07, 0.10]
[0.03, 0.19]
[-0.05, 0.08]
[-0.15, -0.02]
[-0.13, 0.04]
[-0.01, 0.16]
[-0.10, 0.14]
[-0.10, 0.11]
[-0.07, 0.13]
[-0.19, 0.09]
[-0.12, 0.14]
[-0.04, 0.12]
[-0.08, 0.05]
[-0.09, 0.04]
[-0.12, 0.05]
[-0.06, 0.10]
[-0.13, 0.11]
[-0.09, 0.12]
[-0.08, 0.12]
[-0.18, 0.08]
[-0.22, 0.05]
[-0.13, 0.13]
[-0.11, 0.10]
[-0.17, 0.06]
[-0.23, 0.04]
[-0.13, 0.14]
[-0.03, 0.08]
[-0.04, 0.06]
[-0.01, 0.08]
[-0.06, 0.07]
[-0.04, 0.08]
[-0.05, 0.10]
[-0.05, 0.08]
[-0.06, 0.06]
[-0.09, 0.06]
[-0.11, 0.04]
[-0.08, 0.04]
[-0.07, 0.03]
[-0.04, 0.05]
[-0.07, 0.05]
[-0.06, 0.06]
[-0.06, 0.06]
[-0.06, 0.04]
[-0.11, -0.02]
[-0.11, 0.01]
[-0.04, 0.07]
[-0.15, 0.00]
[0.03, 0.16]
[-0.07, 0.05]
[-0.19, -0.02]
[-0.12, 0.04]
[-0.12, 0.15]
[-0.06, 0.17]
[-0.19, 0.05]
[-0.10, 0.19]
[-0.08, 0.19]
Slopes
Group Differences in Slopes
The table below shows the differnces in slopes between the groups. These are basically are critical “significance” tests of the cross level interaction between wave and life event groups.
[-0.06, 0.02]
[-0.05, 0.02]
[-0.02, 0.04]
[-0.06, 0.03]
[-0.05, 0.03]
[-0.06, 0.01]
[-0.03, 0.04]
[-0.04, 0.03]
[-0.01, 0.07]
[-0.05, 0.03]
[-0.04, 0.07]
[-0.01, 0.10]
[-0.05, 0.05]
[-0.03, 0.10]
[-0.05, 0.07]
[-0.03, 0.04]
[0.01, 0.09]
[-0.02, 0.05]
[-0.02, 0.07]
[-0.01, 0.07]
[-0.05, 0.07]
[-0.01, 0.10]
[-0.09, 0.02]
[-0.03, 0.11]
[-0.08, 0.06]
[-0.10, 0.03]
[-0.03, 0.08]
[-0.09, 0.03]
[-0.08, 0.06]
[-0.08, 0.05]
[-0.06, -0.00]
[-0.02, 0.03]
[-0.05, 0.00]
[-0.04, 0.03]
[-0.04, 0.03]
[-0.05, 0.03]
[-0.03, 0.03]
[-0.05, 0.02]
[-0.07, 0.02]
[-0.03, 0.05]
[-0.02, 0.03]
[-0.02, 0.04]
[-0.03, 0.02]
[-0.01, 0.06]
[-0.01, 0.05]
[-0.03, 0.03]
[-0.00, 0.05]
[-0.01, 0.04]
[-0.03, 0.03]
[-0.04, 0.02]
[-0.01, 0.06]
[-0.03, 0.04]
[-0.06, 0.00]
[-0.06, 0.03]
[-0.06, 0.02]
[-0.06, 0.07]
[-0.06, 0.06]
[-0.04, 0.08]
[-0.08, 0.08]
[-0.06, 0.08]
Variances
Were there individual differences in slopes and intercepts? How much of the variance do our models explain. Below, you see a very large table of the level 1 and 2 variances as well as teh correlations among level 2 variances.
Extract Samples
It’s no fun to do Bayesian models without extracting the samples. We’ll use these to look at trace plots and posterior distributions of our models terms. Note that the reason we do this is because we only want to keep samples for some of the parameters. We have several thousand samples for several thousand parameters because the parameters also include individual level effects.
The samples take awhile to extract because we have to load in all the models. So instead, we’ll just load in the samples.
Plot Results
Sample Posterior Distributions
First, let’s plot the posterior distributions. I don’t really want to do this for 70 models separately, so I’ll just show you a sample. You can look at my web app if you want to see all models.
Sample Trace Plot
Let’s also look at a trace plot. This is a trace plot of Markov chain from one multiply imputed Bayseian multilevel growth model. In this sample, the examined trait was Extraversion, and the life event was the a moving out of one’s parents’ home. This is a healthy Markov chain that exhibits both stationarity and well-mixing. The iterations represent only post-warmup samples.
Mean Level Trajectories
Next, we’ll check out mean level trajectories. If the groups are different, we should see differences in change. If not, the trajectories will overlap. I plot these with a 95% interval around them.
Sample Plots
Here’s one sample plot that I use in my talk to get people oriented to the graphs.
Extraversion
Agreeableness
Conscientiousness
Neuroticism
Openness
Individual Level Trajectories
Often, we are also interested in individual differences, so we can look at spaghetti plots of the predicted trajectories for each person for each trait x combination.
Sample Plot
Here’s a sample plot I use in my talk to demonstrate how to understand it.
Extraversion
Agreeableness
Conscientiousness
Neuroticism
Openness
Group Slopes
This is a fun one. For these plots, we’ll plot the estimated slope for each group, the 95% HPDI around that estimate, and the distribution of random slopes around that estimate. “Significant” differences will be flagged with color while the rest will be sad and gray. This code is unwieldy, but the result is pretty.
Sample Plots
Here is a sample plot of group differences I’ll use for my talk.
Extraversion
Agreeableness
Conscientiousness
Neuroticism
Openness
All Slopes